% rsts_STAT.M
% compute posterior mean, sd, and confidence interval (HPD)
% based on Metropolis-Hastings output

% Taeyoung Doh
% created: 02/02/2007
%--------------------------------------------------------------------1.0
% house cleaning
%--------------------------------------------------------------------
clear all;
close all;
diary off;

%--------------------------------------------------------------------
% default parameters
%--------------------------------------------------------------------

% call environments and model selection
rsmodel_spec;
loaddata;

%--------------------------------------------------------------------
% prior specification
%--------------------------------------------------------------------
global prior_;
prior_ = feval(fnames_.prior,msel);


nobs    = size(data,1); 

global datasel psel 
% MHRUN
%  length of mhrun defines the number of columns in graph table
%	odd integer  (e.g., 1) for linear posterior draws
%	even integer (e.g., 2) for quadratic posterior draws
 
mhrun = 1;
mhrunid  = num2str(mhrun','%2.0f');

block1   = 1;		% starting block
blockn   = 100;		% ending block
nsim     = 10000;	% number of simulation draws in each block
nm =1000;
hpdprob  = 0.90;	% probability for highest posterior density

% for data density (modified harmonic mean)
hmax    = 20;

%--------------------------------------------------------------------
% model descriptions and basic calculations
%--------------------------------------------------------------------

% parameter names
pnames = feval(fnames_.pnames,msel);

% reset diary file
diaryfile = [path_.result 'stat_mhrun' mhrunid '.log'];
if exist(diaryfile,'file')
	delete(diaryfile);
end
diary(diaryfile);
diary off;

%--------------------------------------------------------------------
% means and standard deviations
%--------------------------------------------------------------------


% load path
  lpath = [path_.para 'mhrunpm' mhrunid '/' ];  
% initialize output
ndraws     = 0;
drawmean   = 0;
drawsqmean = 0;
sumdrawsq  = 0;
mpos       = zeros(blockn-block1+1,1); 
rho_inf    = zeros(nm,4); 
t=1; 

% loop for block
for jblock=block1:blockn

	% load file
	lfile = [lpath 'mhdraw' num2str(jblock,'%04.0f')];
	load(lfile);

	% use only p.d. covariance matrix results for prior draws
	if mhrun==0
		parasim = parasim(find(retcode==0),:);
	end
     
	% number of simulations in each block
	nsimul = size(parasim,1);

	% collect simulations
	drawblock = parasim;
    drawdim   = size(drawblock,2);
  
	% compute sums of x(i) and x(i)^2 to be used to calculate means and s.d.
	drawmean   = drawmean + sum(drawblock,1);
	drawsqmean = drawsqmean + sum(drawblock.^2,1);
	sumdrawsq  = sumdrawsq + drawblock'*drawblock;
	ndraws     = ndraws + nsimul;
end

% compute means and standard deviations
drawmean   = drawmean/ndraws;
drawsqmean = drawsqmean/ndraws;
drawstdd   = sqrt(drawsqmean - drawmean.^2);
drawsig    = sumdrawsq/ndraws - drawmean'*drawmean;

densfac    = mean(postsim); 
 

% check diagonal elements 
for j=1:size(drawsig,1)
	if drawsig(j,j) < 1E-6
		drawsig(j,j) = 1E-6;
	end
end

drawsiginv = inv(drawsig);
drawsiglndet = log(det(drawsig));

%--------------------------------------------------------------------
% HPD Interval
%--------------------------------------------------------------------
drawci = zeros(2,drawdim);

for j=1:drawdim

	drawcol = [];
	ndraws  = 0;

	% loop for block
	for jblock=block1:blockn

		% load file
		lfile = [lpath 'mhdraw' num2str(jblock,'%04.0f')];
		load(lfile);

		% number of simulations in each block
		[nsimul,npara] = size(parasim);
        drawblock = parasim;  
        
		% Read only the j'th column
		drawcol = [drawcol ; drawblock(:,j)];
		ndraws  = ndraws + nsimul;

	end

	% calculate HPD for j'th column
	drawci(:,j) = hpdint(drawcol,hpdprob);

end



%--------------------------------------------------------------------
% data density:  modified harmonic mean by Geweke (1999)
%--------------------------------------------------------------------
p = (.1:.1:.9)';
pcrit = chi2inv(p,ones(length(p),1)*npara);

ndraws  = 0;
suminvlike = zeros(length(p),1);
laginvlike = zeros(hmax,length(p));
gaminvlike = zeros(hmax,length(p));

paracen = drawmean';   % used for Chib and Jeliazkov (2001) 
weight  = 0;

 ssT     = zeros(nobs,4); 
 rhot    = zeros(nobs,1000);
% loop for block
 for jblock=block1:blockn

	% load file
	 lfile = [lpath 'mhdraw' num2str(jblock,'%04.0f')];
	 load(lfile);

	% number of simulations in each block
	 [nsimul,npara] = size(parasim);				

	 paradev  = parasim-repmat(drawmean,nsimul,1);
	 quadpara = sum((paradev*drawsiginv).*paradev,2);

	% simulation loop
	 for j=1:nsimul
         paraj   = parasim(j,:)';  
        
       lnfpara = - 0.5*npara*log(2*pi) - 0.5*drawsiglndet -0.5*quadpara(j) - log(p);%
 		indpara = (quadpara(j)<pcrit);
	 	invlike = exp(lnfpara - postsim(j) + densfac).*indpara;

	 	laginvlike = [invlike' ; laginvlike(1:hmax-1,:)];
	 	gaminvlike = gaminvlike + laginvlike.*repmat(invlike',hmax,1);
		suminvlike = suminvlike + invlike;
        
          if mprior>0 
        
          k = mod(j,100); 
          
            if k==0
          [sssT,rho_inft]=regime_prob(paraj,data);       
           ssT = ssT + sssT;  
           rho_inf(t,:) = rho_inft';
           rhot(:,t) = sssT*rho_inft;
           t  = t+1; 
          end
          end 
      
	end	% simulation loop

    
	ndraws = ndraws + nsimul;
    
end	% loop for blocks

%ssT = ssT/nm; 

meaninvlike = suminvlike/ndraws;
gaminvlike  = gaminvlike/ndraws - repmat((meaninvlike.^2)',hmax,1);


% standard error
suminvlikeerror = gaminvlike(1,:);
 for k=2:hmax
   suminvlikeerror = suminvlikeerror + 2*gaminvlike(k,:)*(1-(k-1)/hmax);
 end

suminvlikeerror = 100*sqrt(suminvlikeerror/ndraws)./meaninvlike' ;

% marginalized data density
mdd = densfac-log(meaninvlike);
md  = mean(mdd); 
mdd = [p mdd];
 

%--------------------------------------------------------------------
% save and print
%--------------------------------------------------------------------


% print log
fmt='%15.8f ';
disp(' ');
disp(' ');
disp(' ');
diary on;
disp(' ');
disp(sprintf('-------------------------------------------------------------------------'));
disp(         'Parameters      Mean            Stdev                   90% CI ');
for i=1:drawdim 
	disp(sprintf(['%s '  repmat(fmt,1,4) ],pnames(i,:),drawmean(i),drawstdd(i),drawci(2,i),drawci(1,i)));
end
disp(       'Marginal Likelihood : Geweke (1999)'); 
disp(  md); 
disp(sprintf('-------------------------------------------------------------------------'));
diary off;

if mprior==1 
% plot posterior expected probability of regime 

ti = 1953:0.25:2006.75;
plot(ti,ssT(:,1));
title('Posterior Expected Value of Probability of Active Regime'); 
tp1 = corr(data(1:105,2),data(2:106,2));tp2 = corr(data(125:end-1,2),data(126:end,2)); 
mrho_inf = mean(rho_inf)'; stdrho_inf = std(rho_inf)'; 
hpdrho1 = hpdint(rhot',0.9);
figure;
plot(ti,[hpdrho1(2,:)' hpdrho1(1,:)' mean(rhot,2)]); 

title('Time Series Plot of Model Implied Inflation Persistence'); 

elseif mprior==2 

% plot posterior expected probability of regime 

ti = 1953:0.25:2006.75;
figure;
plot(ti,ssT(:,2));
title('Posterior Expected Value of Probability of Low Volatility Regime'); 

tp1 = corr(data(1:105,2),data(2:106,2));tp2 = corr(data(125:end-1,2),data(126:end,2)); 
mrho_inf = mean(rho_inf)'; stdrho_inf = std(rho_inf)'; 
hpdrho1=hpdint(rhot',0.9);
figure; 
 plot(ti,[hpdrho1(2,:)' hpdrho1(1,:)' mean(rhot,2)]); 
 title('Time Series Plot of Model Implied Inflation Persistence');

elseif mprior>=4 
% plot posterior expected probability of regime 
figure; 
ti = 1953:0.25:2006.75;
subplot(2,2,1); plot(ti,ssT(:,1));
subplot(2,2,2); plot(ti,ssT(:,2)); 
subplot(2,2,3); plot(ti,ssT(:,3)); 
subplot(2,2,4); plot(ti,ssT(:,4)); 

 figure; 
 subplot(2,1,1); plot(ti, ((ssT(:,1)+ssT(:,2)))/14); 
 subplot(2,1,2); plot(ti, ((ssT(:,1)+ssT(:,3)))/14); 

 tp1 = corr(data(1:105,2),data(2:106,2));tp2 = corr(data(125:end-1,2),data(126:end,2)); 
 mrho_inf = mean(rho_inf)'; stdrho_inf = std(rho_inf)'; 
 mdd

  hpd1 = hpdint(rho_inf(1:nm,1),0.9);
  hpd2 = hpdint(rho_inf(1:nm,2),0.9); 
  hpd3 = hpdint(rho_inf(1:nm,3),0.9);
  hpd4 = hpdint(rho_inf(1:nm,4),0.9);
   disp(       'HPD Interval for Inflation Persistence'); 
   disp( sprintf(['%s ' repmat(fmt,1,4) ],hpd1(2), hpd1(1), hpd2(2), hpd2(1), hpd3(2), hpd3(1), hpd4(2), hpd4(1) ) ); 
   hpdrho1 = hpdint(rhot',0.9);
   figure; 
   plot(ti,[hpdrho1(2,:)' hpdrho1(1,:)' mean(rhot,2)]); 
    title('Time Series Plot of Model Implied Inflation Persistence');
end 
 
sfile = [path_.para 'stat_mhrunch1' mhrunid '.mat'];
save(sfile,'drawmean','drawstdd','drawsqmean','drawci','drawsig','ssT','hpd1','hpd2','hpd3','hpd4');
